/*!
 * \file rdmaxclique.c
 *
 * \date Sep 24, 2013
 * \author Ricardo C. Corrêa (correa@lia.ufc.br)
 *
 * \brief Russian doll algorithm for maximum clique.
 */

#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <sys/time.h>

#include <bitmap.h>
#include <graph.h>
#include <grinput.h>
#include <heur.h>

/*!
 * \fn int mcr(Graph * g, int * res)
 *
 * \brief Colors input array of vertices in a decreasing order of weights.
 *
 * \param g The input graph.
 * \param res New ordering of the vertices.
 *
 * \return The number of vertices in the input graph.
 */
int wmcr(Graph * g, int * res, const double * weight);
int wmcr(Graph * g, int * res, const double * weight) {
	int j;

	//initializes the vector of degrees
	for (j = 0; j < g->n; j++)
		res[g->n + j] = j;

	//initializes the vector of ex-degrees
	double wdeg[g->n];
	int lastnode = NODEIDX(g->n-1);
	int nodeidx;
	int u;
	for (j = 0; j < g->n; j++) {
		wdeg[j] = 0;
		nodeidx = 0;
		u = NODELEAST(NEIG(g, j)[nodeidx]);
		while (u < 0 && nodeidx < lastnode)
			u = NODELEAST(NEIG(g, j)[++nodeidx]);
		while (u >= 0) {
			wdeg[j] += weight[NELEM(nodeidx) + u];
			u = NODENEXTLEAST(NEIG(g, j)[nodeidx],u);
			while (u < 0 && nodeidx < lastnode)
				u = NODELEAST(NEIG(g, j)[++nodeidx]);
		}
	}

	inline int comparKey1(const void * a, const void * b) {
		return weight[*((int *) b)] < weight[*((int *) a)] ? -1 :
				weight[*((int *) b)] == weight[*((int *) a)] ?
						wdeg[*((int *) b)] < wdeg[*((int *) a)] ? -1 : 1 : 1;
	}

	// sorts vertices according to degrees
	qsort(res + g->n, g->n, sizeof(int), comparKey1);
	for (j = 0; j < g->n; j++)
		res[res[g->n + j]] = j;

	return g->n;
}

typedef struct {
	int vidx;
	int minRNode; //!< a lower bound for the node index of the smallest element in the set
	int maxRNode; //!< an upper bound for the node index of the greatest element in the set
	NODETYPE * R;
	int minWNode; //!< a lower bound for the node index of the smallest element in the set
	int maxWNode; //!< an upper bound for the node index of the greatest element in the set
	NODETYPE * W;
	double maxwcq; //!< Weight of a maximum clique in this doll
} Doll;

static double (*coloring)(NODETYPE *, NODETYPE *, double *, double, int *,
		int *, int *, int *);

void rd_maxDoubleClique(Graph * gg, const double * w, size_t * ndolls,
		size_t * memsz) {
	int n = gg->n;
	int nnodes = NODEIDX(n-1) + 1;

	int * R = calloc(n << 1, sizeof(int));
	wmcr(gg, R, w);
	Graph * g = newGraph(n);
	*memsz = n * sizeof(int);

	Doll * stack = (Doll *) calloc(n, sizeof(Doll));
	Doll * doll = stack;
	double * rw = (double *) calloc(n, sizeof(double));
	*memsz = n * (sizeof(Doll *) + sizeof(double));

	// construction of renumbered graph
	int i, j, k, l;
	for (i = 0; i < n; doll++, i++) {
		doll->R = (NODETYPE *) calloc(nnodes, WSIZE);
		doll->W = (NODETYPE *) calloc(nnodes, WSIZE);
		stack->R[NODEIDX(i)] = ADDTONODE(stack->R[NODEIDX(i)], IDXINNODE(i));

		NODETYPE * ignode = NEIG(g,R[i]);
		NODETYPE * iggnode = NEIG(gg,i);
		rw[R[i]] = w[i];
		k = NODEIDX(i);
		j = NODENEXTLEAST(iggnode[k], IDXINNODE(i));
		while (j >= 0) {
			l = R[NELEM(k) + j];
			ignode[NODEIDX(l)] = ADDTONODE(ignode[NODEIDX(l)],IDXINNODE(l));
			NEIG(g,l)[NODEIDX(R[i])] =
					ADDTONODE(NEIG(g,l)[NODEIDX(R[i])], IDXINNODE(R[i]));
			j = NODENEXTLEAST(iggnode[k], j);
		}
		for (k++; k < nnodes; k++) {
			j = NODELEAST(iggnode[k]);
			while (j >= 0) {
				l = R[NELEM(k) + j];
				ignode[NODEIDX(l)] = ADDTONODE(ignode[NODEIDX(l)],IDXINNODE(l));
				NEIG(g,l)[NODEIDX(R[i])] =
						ADDTONODE(NEIG(g,l)[NODEIDX(R[i])], IDXINNODE(R[i]));
				j = NODENEXTLEAST(iggnode[k], j);
			}
		}
	}
	*memsz += (nnodes << (LOGW + 1)) * n;

	doll = stack;

	doll->minRNode = 0;
	doll->maxRNode = nnodes - 1;
	doll->minWNode = nnodes;
	doll->maxWNode = -1;
	int v = 0;

	initHeur(g);

	//initialize with a maximal clique
	doll->maxwcq = maxClique(doll->R, rw, doll->minRNode, doll->maxRNode);

	//build a coloring of the size of the maximal clique found
	double col = coloring(doll->R, doll->W, rw, doll->maxwcq, &doll->minRNode,
			&doll->maxRNode, &doll->minWNode, &doll->maxWNode);

	printf(
			"Clique weight: %lf  Coloring weight: %lf  minRNode=%d  maxRNode=%d\n",
			doll->maxwcq, col, doll->minRNode, doll->maxRNode);

	if (doll->minRNode > doll->maxRNode)
		goto TERMINATE;

	int * curC = R; //!< Current clique
	int * bestC = R + n; //!<
//	int bottomC = 0;
	int topC = 0;
//	double wC = 0.0;
	int sizeBestC = 0;
	int vnode;
	(*ndolls) = 1;

	while (topC >= 0) {
		//get the next vertex in the candidate set
		doll->vidx = NODELEAST(doll->R[doll->minRNode]);
		if (doll == stack && doll->vidx < 0) {
			doll->minRNode = minNode(doll->R, doll->minRNode, doll->maxRNode);
			if (doll->R[doll->minRNode] == 0)
				goto TERMINATE;
			doll->vidx = NODELEAST(doll->R[doll->minRNode]);
		}

		if (doll->vidx < 0) {
			doll--;
//			if (topC == bottomC)
//				bottomC--;
//			wC -= rw[curC[--topC]];

			if ((doll + 1)->maxwcq + rw[curC[--topC]] > doll->maxwcq) {
				//eliminate a maximal stable set that contains sub->v
//				while (bottomC > 0 && doll->minRNode <= doll->maxRNode)
//					maxStab(doll->R, doll->W, curC[bottomC--], &doll->minRNode, &doll->maxRNode, &doll->minWNode, &doll->maxWNode);
//				if (doll->minRNode <= doll->maxRNode)
				maxStab(doll->R, doll->W, curC[topC], &doll->minRNode,
						&doll->maxRNode, &doll->minWNode, &doll->maxWNode);

				//expand the current best clique
				double expc = (doll + 1)->maxwcq + rw[curC[topC]];
				if ((doll + 1)->maxwcq == 0.0)
					sizeBestC = 0;
				bestC[sizeBestC++] = curC[topC];
				if (doll->minRNode <= doll->maxRNode) {
					expc = expClique(doll->R, rw, bestC, doll->minRNode,
							doll->maxRNode, sizeBestC);
					coloring(doll->R, doll->W, rw, expc - (doll + 1)->maxwcq,
							&doll->minRNode, &doll->maxRNode, &doll->minWNode,
							&doll->maxWNode);
				}

				doll->maxwcq = expc;
			}
		} else {
			v = doll->vidx + NELEM(doll->minRNode);
			if (topC == 0)
				printf("v = %3d(%d) best = %2lf ndolls = %3d\n", v, g->n,
						doll->maxwcq, *ndolls);
			vnode = doll->minRNode;
			doll->R[vnode] = DELFROMNODE(doll->R[vnode], doll->vidx);
			doll->W[vnode] = ADDTONODE(doll->W[vnode], doll->vidx);
			if (vnode > doll->maxWNode)
				doll->maxWNode = vnode;
			if (vnode < doll->minWNode)
				doll->minWNode = vnode;
			doll->minRNode = minNode(doll->R, doll->minRNode, doll->maxRNode);

			//Calculate the new candidate vector
			Doll * subdoll = doll + 1;
			intersectOf(subdoll->R, NEIG(g, v), doll->W, doll->minWNode,
					doll->maxWNode + 1);
			subdoll->minRNode = minNode(subdoll->R, doll->minWNode,
					doll->maxWNode);
			subdoll->maxRNode = maxNode(subdoll->R, doll->maxWNode,
					subdoll->minRNode);

			//Check if the candidate vector is not empty
			if (subdoll->minRNode <= subdoll->maxRNode) { // && cardOfSet(subdoll->R, subdoll->minRNode, subdoll->maxRNode+1) >= max-topC-1) {
				DELALL(subdoll->W+subdoll->minRNode,
						subdoll->maxRNode-subdoll->minRNode+1);
				subdoll->minWNode = nnodes;
				subdoll->maxWNode = -1;

				//Build a coloring for the candidate's vector
				subdoll->maxwcq = doll->maxwcq - rw[v];
				coloring(subdoll->R, subdoll->W, rw, subdoll->maxwcq,
						&subdoll->minRNode, &subdoll->maxRNode,
						&subdoll->minWNode, &subdoll->maxWNode);
			} else
				subdoll->maxwcq = 0.0;
//			if (bottomC == 0 || (topC == bottomC && !doll->R[doll->minRNode]))
//				bottomC++;
			doll = subdoll;
			curC[topC++] = v;
//			wC += rw[v];
			(*ndolls)++;
		}
	}

	TERMINATE:

	gg->maxwcq = stack->maxwcq;

	*memsz += getMemSize();
	for (doll = stack, i = 0; i < n; doll++, i++) {
		free(doll->R);
		free(doll->W);
	}

	closeHeur();
	delGraph(g);
	free(g);
	free(stack);
	free(R);
	free(rw);
}

int main(int argc, char *argv[]) {
	char name[200];
	FILE* graphFile = NULL;
	FILE* weightsFile = NULL;
	FILE* outFile = NULL;

	outFile = fopen(argv[1], "a+");
	graphFile = fopen(argv[2], "r");
	weightsFile = fopen(argv[3], "r");
	strcpy(name, argv[2]);

	Graph * g = readDimacsGraph(graphFile);
	double * w;

	int algo = atoi(argv[4]);
	switch (algo) {
	case 1: {
		coloring = greedyColoring;
		w = readIntWeights(g->n, weightsFile);
		break;
	}
	case 2: {
		printf("Not implemented.\n");
		exit(0);
		break;
	}
	case 3: {
		printf("Not implemented.\n");
		exit(0);
		break;
	}
	case 4: {
		coloring = fracColoring;
		w = readIntWeights(g->n, weightsFile);
		break;
	}
	default: {
		printf("Algorithm choice is invalid.\n");
		exit(0);
	}
	}

	size_t memsz;
	size_t ndolls;

	struct timespec start, end;
	struct timeval tvs, tve;
	float totalTime;

	gettimeofday(&tvs, NULL);
	rd_maxDoubleClique(g, w, &ndolls, &memsz);
	gettimeofday(&tve, NULL);
	start.tv_sec = tvs.tv_sec;
	start.tv_nsec = tvs.tv_usec * 1000;
	end.tv_sec = tve.tv_sec;
	end.tv_nsec = tve.tv_usec * 1000;
	totalTime = (end.tv_sec - start.tv_sec);
	totalTime += (end.tv_nsec - start.tv_nsec) / 1000000000.0;

	printf(
			"%s \t n=%d  option=%d  sol=%lf  ndolls=%d  time=%f  mem=%fKB  graphmem=%fKB \n",
			name, g->n, algo, g->maxwcq, ndolls, totalTime, memsz / 1024.0,
			g->matrixsize / 1024.0);
	fprintf(outFile, "%s \t %d \t %d \t %lf \t %d \t %f \t %f \t %f\n", name,
			g->n, algo, g->maxwcq, ndolls, totalTime, memsz / 1024.0,
			g->matrixsize / 1024.0);
	fclose(outFile);
	fclose(graphFile);
	free(w);
	delGraph(g);
	free(g);
	return 0;
}
